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Abstract 

At the Haverah Park Array a number of air shower observables were measured that are relevant to 
^ I the determination of the mass composition of cosmic rays. In this paper we discuss measurements of 

the risetime of signals in large area water-Cherenkov detectors and of the lateral distribution function 
of the water-Cherenkov signal. The former are used to demonstrate that the CORSIKA code, using 
the QGSJET98 model, gives an adequate description of the data with a low sensitivity, in this energy 
■ range, to assumptions about primary mass. By contrast the lateral distribution is sufficiently well 

' measured that there is mass sensitivity. We argue that in the range 0.2-1.0 EeV the data are well 

, represented with a bi-modal composition of (34 ± 2)% protons and the rest iron. We also discuss the 

' systematic errors induced by the choice of hadronic model. 

O ■ 1 Introduction 
^ . . 

C/3 . High-energy cosmic rays are measured via extensive air showers of secondary particles they produce in the 

Earth's atmosphere. Specific observables depend in a complex way on primary mass and energy, which 
must be understood to interpret air shower data. Some measurable parameters at ground can be used to 
obtain the primary energy of the incoming cosmic rays without large systematic uncertainties due to the 
unknown composition ||l[ . The main effect of a variable mass on showers with the same initial energy, 
I E, is a characteristic shift of the position of the maximum of shower development in the atmosphere, 

X^. Cosmic ray nuclei share their energy amongst A nucleons and, therefore, their showers can, to a 
first approximation, be described as a superposition of A nucleon-induced subshowers with an energy of 
E/A each. Thus, showers of heavy primaries are less penetrating, and tend to develop higher up in the 
atmosphere, than nucleon showers of the same energy. Therefore, is an observable with a strong 
connection to the primary mass. However, additionally, the depth of maximum also depends on energy. 
The higher the energy the longer the shower and the deeper is X^. This is described by the elongation 
rate, dX^/d\ogE, which is usually quoted as the shift in Xm per change of energy by a factor 10. In 
practice, apart from the unique case of fluorescence detectors, there is no way to measure X^ directly. 
Mostly other quantities are recorded that can be shown to have a correlation to the height of the shower 
development. 

Here we report an analysis of mass composition in the energy range 0.2-1.0 EeV that has been performed 
with data from the Haverah Park extensive air shower array. The Haverah Park array was a 12 km^ air 
shower array consisting of water tanks that acted as Cherenkov detectors. It was operational from 1967 
to 1987. 

At Haverah Park a number of observables were measured which are relevant to the determination of the 
mass composition of cosmic rays above 2 x 10^^ eV. In particular, the two parameters, rj and tj^/2: which 
are sensitive to the longitudinal development of showers, have been studied in detail, rj describes the 

*now at: Max-Planck-Institut fiir Kernphysik, Heidelberg, Germany 



1 







1 00 rr 




A2 
• 






o 




o 

o 

o 

□ 

A 4 

4 A1 o • 
o • * □ 


o 


1 m 
2.3 




□ 


□ 
• 


9 

o 

34 m' 


o 

o 

o 




A3 
• 


o 



Figure 1: The inner part of the Haverah Park array, the so called infill array. 

steepness of the lateral distribution function of the signal observed in the water-Cherenkov detectors, rj 
was measured with high precision using a portion of the array with a much denser detector arrangement, 
the so-called infill array For 1425 events in the energy range 2 x 10^^ to 3 x 10^^ eV information from 
infill detectors was available and 77 could be determined. These were recorded between 1977 and 1981. 
The second mass-sensitive parameter, the risetime ti/2, is a measure characterising the spread of the 
arrival times of individual particles at a given detector. It is evaluated from the integrated signal of all 
particles belonging to a shower and defines the time interval in which the integrated signal rises from 10% 
to 50%. Risetimes can be sensibly determined only in large-area detectors. Therefore they are obtained 
from the four 34 m^ detectors, A1-A4 (Fig. |^). Risetime data were obtained from > 7000 events and a 
total of 13000 detector signals at core distances of more than 300 m. The events used here span a range 
of energies from 2 x 10^'' to 2 x 10^^ eV. 

The analysis of risetimes of events with known core position, arrival direction and energy established 
the first evidence of shower-to-shower fluctuations at high energies 

[§• This was supported by the 

demonstration of a correlation of ti/2 with the steepness of the lateral distribution, which could be 
understood as a consequence of both variables being dependent on X^^- In the 1970s and 1980s it was 
not possible to make accurate predictions of the ti/2 expected for different primary masses because the 
problem could not be solved analytically and a full 4-dimensional numerical simulation was beyond the 
capabilities of the computing power accessible at that time. 

The analysis of 77 data also showed evidence of fluctuations very much larger than could be attributed 
to the experimental uncertainties The dependence of rj on the depth of the shower development 
in the atmosphere could not be established directly from the data. Furthermore, a comparison of the 
measured average value of 77 with a highly regarded model calculation of the time ||] showed strikingly 
poor agreement: a mean cosmic-ray mass very much heavier than iron was required to fit the data (see 
Fig. |). 

The situation is now very different and a model can be found that describes the data rather well. To 
interpret the data detailed simulations with the CORSIKA and GEANT programs have been performed. 
CORSIKA 1^ is a modern air shower model that simulates particle production and propagation in the 
atmosphere fully. It employs various models of high-energy hadronic interactions, of which the most 
tested is QGSJET98 (Quark-Gluon-String model with jets) 0. The GEANT program § is used to 
account for the detector response to particles of different type, energy and impact angle. 
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Figure 2: A comparison of Haverah Park data and the calculations for p and Fe induced showers from |^]. p is 
the water-Cherenkov signal. 

2 The Haverah Park Array 

The Haverah Park array has been described in detail elsewhere The central part of the array is 
composed of four 34 water-Cherenkov detectors, A1-A4, spaced 500 m apart. Six groups of satellite 
arrays each of 4 x 14 lie near to a circumference of radius 2 km around the centre so that density 
information is available from all sites for showers that trigger the 500 m array. The 500 m array is 
sensitive to showers with primary energies above w 6 x 10^^ eV. Another set of three detectors (9 
each) surrounding the central detector (Al) form a 150 m array. The so-called infill array consists of 23 
small-area water-Cherenkov detectors located in the central region of the 500 m array with an irregular 
spacing of « 150 m. For the arrangement of detectors, see Fig. 0. Further small-area detectors were 
placed right next to the 9 w? detectors and Al for calibration purposes. 

There were two kind of detectors: (i) 2.29 vc? area galvanized rectangular iron tanks with aluminum lids, 
identical to those used in the 4 x 34 detectors: these detectors were placed inside wooden huts with 
light roofs, and (ii) 1 vc? area water tanks, made of light-weight polyurethane foam and covered with fibre 
glass: these tanks of hexagonal cross section, stood outdoors. All tanks were filled with water to a depth 
of 1.2 m and viewed by one photomultiplier with 100 cm^ photocathode. Detector areas larger than 2.29 
m^ were achieved by grouping several of the 2.29 rv? detectors together. Sixteen of the modules form 
each of the detectors A1-A4. The signals from 15 of the 16 tanks were summed to provide the signal 
used for triggering and for the density estimate. The 16"^ tank in each group was used to provide a low 
gain signal. 4 modules make up the 9 w? detectors and 6 modules form the 13 ir? detectors. 
The densities of Cherenkov photons per unit detector area (Cherenkov densities) were recorded in terms 
of the average signal from a vertical muon (1 vertical equivalent muon — 1 vem) per square metre. This 
signal corresponds to approximately 14 photoelectrons (pe) for Haverah Park tanks. A trigger for all 
detectors is generated when the central detector (Al) and at least two out of the three other A-sites (A2, 
A3, A4) detect a particle density of > 0.3 vem/m^. However, infill signals had to be above 7 vem/m^ 
to be recorded. Trigger rates were monitored daily over the life of the experiment. After correction for 
atmospheric pressure variations, the trigger rates were found to be stable to better than 5%. 
A feature of the Haverah Park array, which is particularly important for the shower front measurements, 
is the area of the detectors A1-A4. These detectors allow large samples of the shower front to be taken 
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at widely separated points, thus minimizing the chance of observing effects due only to local density 
fluctuations. The temporal response of the recording system to a step function was determined by 
measuring the response to small showers that fell nearby and thus produced short risetimes. The risetime 
of bandwidth limited pulses was 45 ns 0. 

Cross-calibration of the detectors of the infilled array with the other detectors could have been achieved 
using vertical muons, as for the large-area detectors, so that they yield the same response to a purely 
muonic signal. However with such a cross-calibration one would have had to correct for the different 
response of the detectors to the soft component and for differences in optical and geometrical properties 
of the detectors and so an alternative approach was adopted. Signals from individual 2.29 m^ and 1 
m^ detectors, next to Al, were compared with corresponding densities from the 34 m^ detector, and an 
appropriate conversion factor was obtained. Similar calibration procedures were used on simulated data 
for consistency. 



3 Model calculations 

At present the CORSIKA program is effectively the standard tool for air shower simulations, and it is 
successfully used by a variety of different experiments over the energy range from 1-10^^ GeV. CORSIKA 
uses EGS4 (Electron Gamma Shower code) for the simulation of the electromagnetic interactions, 
and features a detailed simulation of particle propagation and decay. However, most important for 
the development of air showers in the atmosphere are soft hadronic interactions. At present the only 
available theoretical approach to model these is Gribov-Regge theory of multi-Pomeron exchange. Several 
interaction models based on this theory are available in CORSIKA, the most successful being QGSJET 
[0. The QGSJET model also contains the treatment of hard collisions and the production of mini-jets, 
which become dominant at very high energy collisions, and includes a realistic simulations of nucleus- 
nucleus collisions. Its parameters were tuned to fit a wide range of accelerator results and to provide 
a secure extrapolation to the highest energies. The combination of CORSIKA with QGSJET seems 
to be able to describe experimental cosmic ray results from Cherenkov telescopes at 10^^ eV over 
measurements of the hadronic, muonic and electromagnetic shower components in the knee region |13|, 
to lateral distributions and arrival times as will be shown in this work, up to air showers at 10^° eV |14|. 
Using the CORSIKA code with QGS JET98 we have generated a library of proton and iron showers with 
zenith angles: 0°, 15°, 26°, 40° and 45° at a primary energy of 4 x 10^^ eV, and for a zenith angle of 
26° at energies of 0.2, 0.4, 0.8, 1.6, and 3.2 EeV. For each set 100 showers with the same parameters have 
been simulated. Additionally we have generated sets of 100 showers with oxygen and helium primaries at 
zenith angle 26° and energy 0.4 EeV. In total 3600 showers were produced. Statistical thinning of shower 
particles was applied at the level of 10~^x E with a maximum particle weight limit of 10"^'^ x E/eV. 
Predictions of the risetimes were deduced from the time information of each individual particle contained 
in the simulated shower library. The time distribution of the signals produced by electromagnetic particles 
for different core distances (r) was obtained for each simulated shower. The integral over the distribution 
was normalized to the expected signal from electromagnetic particles at the given distance for a 34 m^ 
detector. The arrival time distribution of muons has been obtained for different r. The expected number 
of muons at a given r for a 34 m^ detector has been computed and Poisson fluctuations added. The 
muons are sampled from the arrival time distribution to obtain the time distribution of the muon signal 
in the detector. The time distributions of the muon signal and the soft component were then added and 
convoluted with the known system response to an instantaneous pulse. The risetime is then inferred from 
the result of the convolution. This procedure was repeated 100 times for each shower to get the mean 
risetime as a function of r and the expected spread of the risetimes arising from Poisson fluctuations in 
the number of muons. To increase the statistics further each CORSIKA shower was thrown 100 times 
onto the array with random core positions ranging out to 300 m from the centre of the array. The risetime 
at each of the four 34 m^ detectors is calculated from the parameterisations obtained above. The risetime 
at each detector is smeared according to an error that contains contributions from Poisson fluctuations 
in the number of muons and from the experimental measurement error (« 3 ns). The core distance of 
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each detector was varied according to the error in the core position described in |^]. To reproduce the 
multiphcity of measured risetimes per event in the data, we keep only the core position and the risetime 
of some of the detectors. We repeated this procedure for each of the 100 CORSIKA showers in a set and 
the mean risetime versus r is obtained (Fig. ^. 

To obtain a value of r] for each simulated shower we have convoluted the CORSIKA output with the 
detector response and fitted the resulting lateral distribution function for muons and electromagnetic 
particles. Again, each shower was used 100 times with core positions randomly scattered over the infilled 
area, obtaining the densities at each detector with the parameterisations described above. The densities 
are modified according to Poisson fiuctuations in the number of particles. We tested whether an event 
meets the array trigger conditions and, if so, the densities were fluctuated according to measurement 
errors and recorded in the same format as real data. The calibration of the infilled array is reproduced 
with simulations to obtain the conversion factor described in the previous section. To obtain the energy of 
the showers, either simulated or real, we need its relation to the water-Cherenkov signal density at a core 
distance of 600 m, p(600), and the attenuation length A. The procedure to calculate these is described 
in ref. ||l|. Here we have followed the same procedure. The uncertainty in the core position, for the rj 
analysis, is « 5 m for all energies and primary masses tested. The energy resolution for proton primaries 
falls from 15% at 0.4 EeV to 10% at 6.4 EeV, while for iron primaries it is 10% and 7%, respectively. 
These uncertainties include physical fluctuations in p(600) and measurement errors. 



4 Risetime analysis 

The dependence of the risetime of the signal in a detector on shower development arises largely because 
of geometry. Particle scattering, velocity differences and geomagnetic deflections are second order effects: 
the geometrical sensitivity is enhanced at large distances. The set of nuclear interactions that produce 
particles observed at an off-axis detector may be regarded as lying on a line source. There are two 
important properties of this line source that relate to mass composition and affect the observed risetime. 
These are its position in the atmosphere and its length. 

Previous studies of ti/2 have revealed its dependence on 6, energy and core distance, and used deviations 
from mean values as measures of the variations of with energy to get results on mass composition. 
The variations deduced [p^ were in good accord with measurements made subsequently and directly by 
the Fly's Eye group [l^. 
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Figure 3: Risetime versus distance to shower core as obtained in data and in simulations for proton and iron 
primaries. 
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The strategy for the present analysis is different: the risetime as a function of core distance for different 
bins of 6 and E are compared directly with predictions from Monte Carlo simulations for different primary 
masses. In the energy and distance ranges (250 < r < 500 m) of interest here, where 77 has been obtained 
with high accuracy, the sensitivity of the risetime technique for the extraction of mass information is 
rather limited. The technique is, however, expected to be very effective for mass separation at greater 
distances and in larger showers where the mean risetime difference predicted between proton and iron 
showers is large and easy to measure. The results of such a study will be described elsewhere. At small 
distances limitations of bandwidth mean that relatively large samples of showers were required to establish 
the dependence of risetime on energy and zenith angle and to demonstrate shower-to-shower fluctuations 
1^ , ^ . However this relative insensitivity to mass can be turned to advantage by using the data to test 
whether a shower model is able to predict average risetimes accurately over the distance range discussed 
here. In Fig. |^ we compare the mean risetime with distance, for two zenith angle bands, with model 
predictions. The data lie between the predictions for protons and iron: it is not wise to deduce any mass 
information from these plots because of known systematic uncertainties in the data and the models at 
the few nanosecond level. The flattening-off evident at the smaller core distances is a consequence of the 
low bandwidth of the 1970s recording system. 

Two cuts were applied to the data and to the simulation results to obtain Fig. |^: the density in a 
particular detector had to be greater that 1 vem/m^ (to reduce sampling fluctuations), and the position 
of the core had to be closer than 300 m from the central triggering detector (to avoid large core errors). 
We infer from this analysis that the CORSIKA code, using QGSJET98 physics, gives a very adequate 
description of the mean shower risetime as a function of distance from 200-800 m and over the zenith angle 
range 15°-40°. This is the first time that such agreement has been demonstrated. While the comparisons 
of Fig. H do not provide proof that the QGS JET98 model is correct they do give us reasonable confidence 
in using this model to attempt to interpret other data that are expected to show mass sensitivity. 



5 Lateral distribution analysis 
5.1 Event reconstruction 

First the shower direction was determined using the arrival time information of the four central triggering 
detectors. The particle density information was then analysed to find the shower core. The lateral 
distribution, parameterised with respect to the known core position and the primary energy, is estimated 
from the particle densities and the form of the lateral distribution using the correlation between p(600) 
and energy. At a given energy the lateral distribution was found to be well described experimentally by 
the modified power law: 

p{r) = k r-(''+'-/«"0 -) (1) 

and fits well to the average data in the distance range 80-800 m. r/ was shown to vary with zenith angle 
as 77 3.78 - 1.44 sec 6*. 

Due to the small number of densities usually available without the infill array, it was not possible to 
fit reliably values of t] for each individual shower. Therefore, the dependence of 77 on zenith angle was 
obtained by calculating the average water-Cherenkov lateral distribution function for different zenith 
angle bins. However, the great number of densities available with the infill array allows study of the 
fluctuations of the lateral distribution on a shower-by-shower basis. 

The original algorithm used to determine the shower parameters is described as follows. For a given 
lateral distribution function geometrical reconstruction provides a unique core location using only a 
small number of density measurements. For each selected event the core was found using 3 or 4 detectors 
surrounding the probable core position. The geometrical mean distance of the shower core from these 
"ringing" detectors was in the range 50 < rn < 200 m. Only densities at distances > were used to 
determine the lateral distribution function using a minimisation technique. An iterative procedure 
was used to find the best value of 77 for each event using as a starting point the energy independent 
approximation of 77 given above. A multiple linear regression was used to determine the dependence of 77 
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Figure 4: Example of the reconstruction of an infill event. Left panel: Projection of the array into the shower 
plane with recorded densities shown as circles with radius proportional to the logarithm of the density. The 
detector areas are indicated by grey scales. The detectors are displayed in the plane perpendicular to the shower 
axis. Right panel: Fitted lateral distribution function. The two lines correspond to the lateral distribution 
function as obtained with ringing analysis and with the method described in the text. The difference in the value 
of ri for this particular event is 0.04. The abscissa shows a quantity that linearises the lateral distribution function 
given in Eq. 0. Filled symbols indicate detectors with signals above saturation or below threshold. 

on zenith angle and energy. The resuh was 

r] = a-b (sec 61-1) + clog(£;/10"eV) 

with a = 2.20 ± 0.01, 6 = 1.29 ± 0.05, and c = 0.165 ± 0.022. The difficulty of the ringing analysis 
is the implementation of the algorithm. The optimum set of ringing detectors was chosen subjectively. 
Several factors were considered in the choice of ringing analysis but the experience of the observer was 
also important. For the purposes of this analysis, in which a large statistical sample of simulated events 
is generated to compare with data, another technique described below was adopted to find the optimum 
shower parameters, so that the role of the observer was removed. 

The algorithm used here is a grid search over many different core position. For a given core position 
the best value of 77 and k are found through an analytical fit and the function is computed. We used 
only detector densities above threshold and below saturation in the range 80-800 m. Densities above 
saturation and below threshold do not improve the fit because of the large number of recorded densities 
available in this data set. Fig. ^ shows an example of a reconstructed event using this method and the 
ringing analysis. The difference in the value of 77 between the two analyses were found to be ~ 0.08, 
which agrees with the measurement error reconstruction in 77 for the two techniques (~ 0.08). 
The selection criteria applied to the data are the following: 

• Zenith angles in the range 0°-45°. 

• In every shower the core must be located within the infilled area. 

• Showers in which the largest density is at one of the boundary detectors are excluded. 

• The value of from the 77 fitting procedure should correspond to the "goodness of the fit" having 
a probability greater than 1%. 
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After the application of these criteria 1351 showers remain from the initial number of 1450. This is the 
data set we use to compare with model calculations. 



5.2 Model comparison 

We use two independent methods to extract information about mass composition from the lateral distri- 
bution function. The first one is based on the method and the second one on the likelihood approach. 
The first technique is more direct and allows a visual check of the procedure, but requires the use of 
zenith angle cuts in our data set, reducing thus the statistics. The second method allows us to combine 
all the zenith angles to predict the mass composition in different energy ranges. We will show that the 
two methods are consistent but we will give our final results in terms of the second method. 



5.2.1 Method 1 

The simulated events generated from CORSIKA showers are analysed with the same algorithms as real 
data and the same cuts are applied. In Fig. ^(a) the variation of ij with zenith angle is shown for events 
with energies in the range 0.3-0.5 EeV. It is evident that, if the mean mass lies between the limits of 
proton and iron, the lower energy data are well described by the QGSJET98 model. 
In Fig. ||(b) and ^ we show the variation of rj with energy. In Fig. ^(b) data in three zenith angle ranges 
are displayed. The model calculations shown for the four mass species are for the boundary between the 
two most vertical zenith angle bands. In Fig. |^ the data have been normalised to 26° for comparison 
with calculations at the same angle. The agreement between the normalised data and a smaller data set 
in the range 1.06 < sec 9 < 1.16 is good and gives confidence in the normalisation procedure. A linear fit 
to the normalised data is shown: the reduced is 1-1. The data might be described with a single mass 
component, independent of energy, or by a mixture of several mass components. To distinguish between 
these possibilities we use data on the spread of ij. 

Although the average properties of showers do not vary strongly with primary mass, fluctuations in the 
longitudinal development are substantially larger for proton than for iron-initiated showers. The smaller 
fluctuations for nucleus-induced showers is a direct consequence of the superposition of A nucleonic 
showers, produced after the initial nucleus fragments, which tends to average out extreme fluctuations that 
might occur in individual subshowers. To study fluctuations in 77, two rj distribution were constructed: 
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Figure 5: Right panel: evolution of rj with zenith angle; the simulation results correspond to a zenith angle of 26°, 
the data is binned in zenith angle bands. Left panel: evolution of rj with zenith angle for data and simulations. 
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Figure 6; Evolution of 77 with energy for data and simulation results. 
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(i) showers with 1.06< sec 6* < 1.16 and 0.2 EeV< E < 0.6 EeV (292 events), and (ii) showers with 
1.04 < sec6l < 1.18 and 0.6 EeV< E < 1.0 EeV (46 events). 

To compare with simulations, we account for aU sources that could contribute to the spread of 77, namely 
fluctuations in the longitudinal development of the shower and reconstruction errors. These are automat- 
ically taken into account, since fluctuations are contained in the simulated event sample, and the same 
reconstruction methods are used for simulated events as for real ones. The typically root mean square 
spread from "shower-to-shower fluctuations" is 0.12 for proton showers and 0.05 for iron showers, and 
from the measurement error in r] is « 0.08. 

In Fig. 0(a) and (b) the experimental data in the energy band 0.2-0.6 EeV are compared with predictions 
for proton and iron beams. The reduced for the fits are 31 and 6.4, respectively. The tail at large 77 in 
the comparison with iron indicates that some light nuclei are required to fit the measurements. For helium 
and oxygen the corresponding values are 20.0 and 7.0. A fit to a dual component mixture is shown 
in Fig. 0(c). A mixture with (29 ± 4)% of protons and a corresponding amount of iron is a very good 
fit (x^/dof ~ 1.04): the addition of small amounts of hehum and oxygen (< 2%) gives x^/dof = 0.8. In 
Fig. 0(d) we show 46 events in the energy range 0.6-1.0 EeV and zenith angle range 1.04 < sec^? < 1.18: 
a two component fit gives (30 ± 10)% protons. 

We thus conclude that there is no evidence for any change of mass with energy in the range studied 
and that Fe is the dominant component. The data of Fig. do not exclude the possibility of the Fe 
component getting larger as the energy reaches 1 EeV and then smaller beyond. 



5.2.2 Method 2 

The zenith angle cuts used in the procedure described above do not allow us to use all the data available 
in a given energy range to extract mass composition. For this reason we have developed an additional 
technique, using the likelihood method, that exploits the data more effectively. 

For each set of CORSIKA showers at a given energy and zenith angle we have generated an 77 distribution. 
Each set consists of 100 CORSIKA showers, and each shower was thrown randomly onto the array. We 
have reconstructed it with the same procedure as for the data, obtaining an i] distribution that includes 
the effect of shower-to-shower fluctuations and experimental reconstruction. We have fitted each resulting 
rj distribution to a Gaussian obtaining the mean /i,, and width cr,,. 

We can parametcrisc the values of and cr,, as a function of zenith angle at a fixed energy (0.4 EeV), 
and as a function of the energy at a fixed zenith angle (26°), for proton and iron primaries. Fig. and 
show these fits, which are very satisfactory. Therefore we can predict the values of these parameters for 
any energy and zenith angle. 
The values of /i^ and cr^ are given by: 

firj{0,E) ^ ai + fla Ae + aa A^ + a4 logio(£;/0.4EeV) (2) 
a^{e,E,p) - bi + b2 Ag + bs \og,oiE/OAEeY) (3) 
ar,{0,E,Fe) = b^ + b2 Ag + b3 eMh ^ogwiE/OAEeV)) (4) 

with Ag = sec 9 — sec 26°, and the constants obtained from the fits shown in Figs. and for proton and 
iron primaries are given in Table 0. It should be noticed that we use a different functional form for cr,, 
depending on the primary, cr,, is decreasing with increasing energy as the shower-to-shower fluctuations 
are reduced, but it cannot have a value smaller than the reconstruction error of rj (« 0.08), which is 
constant with energy. This is the reason for the flattening of the value of ct,; for iron primaries at high 
energies in Fig. 0. 

The probability density of an event to belong to a proton or iron i] distribution is given by the appropriate 
Gaussians: 

G = _i exp(-fa-^/'-^;»') (5) 

where 9i, Ei and rji are the zenith angle, the energy and the reconstructed values of rj of each event. 
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Table 1: Parameters obtained in the fits of Eqs. H-ll to simulation results 
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Figure 8: Variation of and ct^ with zenith angle for proton and iron primaries at an energy of 0.4 EeV. The 
lines have the form of Eqs. and are fitted to the results obtained from simulations. The top lines correspond 
to proton primaries and the bottom lines to iron primaries. 
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Figure 9: Variation of and an with primary energy for proton and iron primaries at a zenith angle of 26° 
The lines have the form of Eqs. 0-14, and are fitted to the results obtained from simulations. 



11 



1 

0.9 
0,8 
0.7 
0.6 

0.4 
0.3 
0.2 
0.1 




x'/clof= 10.5/6 



1.05 



1.1 



1.15 
see© 



1.2 



1.25 



1.3 



Figure 10: Predicted value of Fp in the energy range 0.3-0.5 EeV for different zenith angle bins. The predicted 
value of Fp does not depend on zenith angle. 



To fit the composition, assuming only two components as seems reasonable from the discussion in 5.2.1 
we use the the maximum likelihood method. The quantity to maximise in this method is: 



hi PiFp) = ln(Pi P2 .... Pn) - J2 



(6) 



where Pi is given by: 



P 



FpGp 



(1 - Fp)G- 



Fc 



(7) 



where Fp is the fraction of protons and Gp (Gpc) is the probability of the event to be proton (iron), as 
given in Eq. ^. 

With this procedure we can use events at all zenith angles to fit Fp in a given energy range. As a 
demonstration of reliability of the method, we have divided our set of data into bins of zenith angle, 
selecting events between 0.3-0.5 EeV. For each angular bin we fitted the composition using the method 
described above. Fig. |o| shows the results. As expected, the values derived of Fp obtained do not 
depend on zenith angle. The average value is in agreement with the value obtained in the previous 
section (Fp = 0.30 ±0.04). 

We have also investigated the possibility of having oxygen and helium in the cosmic ray beam. For this 
purpose we have selected events in the energy range 0.3-0.5 EeV and with 1.1 < sec 9 < 1.14 , and, using 
a modification of Eq. ^ fitted the fraction of oxygen (helium) assuming a three component mixture of 
primaries: proton, iron, and oxygen (helium). The fitted fraction of oxygen or helium are less than 1%. 
At 68% confidence level the fraction of oxygen (helium) should be less than 15% (25%). Therefore, as in 
the previous section, we will assume a dual component mixture of proton and iron for our final discussion. 
Fig. 11 shows the predicted fraction of protons in the cosmic ray beam as a function of the energy. In 



this case we have used the events from all zenith angles. It is apparent that in the energy range 0.2-1.0 
EeV the composition does not change within the limits of measurement. The rapid increase of Fp at 
energies below 0.2 EeV is an artificial effect. It was established very early (l8| that below 0.2 EeV steep 
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Figure 11: Predicted value of Fp as a function of the energy. A fit to a constant composition in the energy range 
0.2-1.0 EeV is also shown with its corresponding x^- The number of events in each energy bin is shown. The 
shadow region corresponds to the energy range in which the analysis is affected by trigger biases (see text). 




Figure 12: Likelihood curves in two energy bands of Fig. |ll|. The horizontal lines correspond to a value of the 
maximum of the log-likelihood function minus 0.5 (1 a level). 
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showers have a higher probabihty to trigger the array, and as the average proton shower is steeper than 
the average iron shower, Fp increases rapidly at lower energies. A fit to a constant composition in the 
energy range 0.2-1.0 EeV gives a predicted value of = 0.34 ± 0.02. Examples of the likelihood curves 
used to extract mass composition are shown in Fig. |l^ for two energy bins of Fig. |Tl|. 
We thus conclude that there is no evidence for any change of mass with energy in the range studied 
and that Fe is the dominant component in this energy range. The data of Fig. |ll| do not exclude the 
possibility of the proton component being enhanced at energies larger than 1.0 EeV. 

5.3 Sensitivity to hadronic model 

The simulations presented in this work were obtained using CORSIKA and QGSJET98 as high energy 
interaction model. Recently QGSJET98 has been replaced by QGSJETOl. The major change was a 
modification in the treatment of diffractive processes . This change reduced the value of the position 
of the shower maximum by « 10 g/cm^, that is showers develop higher in the atmosphere. 
To establish the sensitivity of our results to this change in the hadronic model we simulated 4 sets of 100 
CORSIKA showers with the new model for proton and iron primaries at a zenith angle of 26° and with 
primary energies of 0.4, and 0.8 EeV. We found that the spread of 77 (cr^) does not change, but the value 
/ir; is reduced by 0.04 for both proton and iron. Therefore, we can calculate the predicted value of Fp as 
a function of the energy using the likelihood method, but reducing the value of /i^ in Eq. ^ by 0.04. Fig. 
|l3| shows the results compared with the expectation obtained in the previous subsection. The value of 
Fp is increased from 0.34 to 0.48. 

The QGSJETOl model has been shown to have the smallest value of when compared to a variety of 
hadronic models (neXus 2, SIBYLL 2.1, DPMJET II. 5) jl^. Some of these models have the mean value 
of X„i increased by ^ 10 g/cm^ compared to QGSJET98. Assuming that the ratio of electromagnetic 
to muonic particles in the shower for these models is the same than QGSJET98, we can estimate the 
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Figure 13: Predicted value of Fp as a function of the energy for QGSJETOl (solid line) and QGSJET98 (dashed 
line). The shadow region corresponds to the energy range in which the analysis is aflfected by trigger biases (see 
text). 
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changes in our results in the manner used for the QGSJETOl model. The value is now increased by 
0.04 and the value of Fp is reduced from 0.34 (QGSJET98) to 0.23. 

We can also estimate the changes in the risetime analysis when using another hadronic model. Using the 
relationship between Xm and ti/2 in a shift of Xm by « 10 g/cm^ corresponds to a change in ti/2 
at 400 m from the shower core of less than 1 ns. Hence, the agreement between data and simulations 
exemplified by Fig. || still holds. 

As a final result, we state that the data on lateral distributions in the energy range 0.2-1.0 EeV is well 
fitted by a dual p/Fe composition with about (34 ± 2)% of the signal being protons constant in this 
energy range. The error in this quantity comes from the statistics of the data available. The systematic 
uncertainty, introduced by the choice of hadronic models, is ~ 14%, larger than the statistical error. 
Clearly, further calculations, using a range of models, are highly desirable. 

6 Comparison with other experiments and conclusions 

The Haverah Park data, in the light of new hadronic models have provided a new estimation on the 
mass composition in the energy range 0.2-1.0 EeV. Efforts to understand the origin of cosmic rays at any 
energy are greatly hampered by our lack of knowledge of the mass distribution in the incoming cosmic 
ray beam. We have obtained here a prediction for the mass composition in the energy range 0.2-1.0 EeV. 
We can compare only with one experiment operating in the same energy range: the HiRes prototype, 
operated with the MIA detector. Using the QGSJET98 model they find a rapid change towards a light 
composition between 0.1-1.0 EeV |2^, which is in contradiction with our results. The question of mass 
composition is thus far from being resolved and there is clearly scope for alternative approaches. 
Some years ago an analysis of the Fly's Eye data j2^ pointed to a change from an iron dominated 
composition at 3 x 10^^ eV to a proton-dominated composition near 10^^ eV. This conclusion was drawn 
from a study of the variation of depth of maximum with energy (the elongation rate) and from an analysis 
of the spread in depth of maximum at a given energy. The conclusions on mass composition from the 
study of the spread in depth of maximum were confirmed by but none of these analyses extended 
to energies above 10"'^^ eV. A different conclusion has been reached by the AGASA group based on the 
variation of the muon content of showers with energy [p^ . Their analysis favours a composition that 
remains "mixed" over the 10^^ to 10^^ eV decade. The statistics in this energy range from this data set 
are too small to be able to claim any agreement with either of the two experiments. We are starting to 
work on a risetime analysis in this energy range which could provide some relevant information. 
We can also compare our results with the mass composition obtained by KASCADE. The mean value of 
the mass at their highest accessible energy (« 10^^ eV) was estimated to be (In A) = 3.5 ±0.5 which 
should be compared with the result obtained in this work which corresponds to (In A) = 2.65 ±0.08 ±0.6 
in the energy range 2-10 x 10^^ eV. In contrast to our results KASCADE also claims that protons have 
largely disappeared at 10^^ eV. 
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